module icosahedron_mod

  use container
  use flogger
  use string
  use const_mod
  use sphere_geometry_mod

  implicit none

  type icosahedron_type
    integer :: nsub = 0 ! Number of subdivision
    integer :: nv   = 0 ! Number of vertices
    integer :: nf   = 0 ! Number of faces or cells
    integer :: ne   = 0 ! Number of edges
    integer, allocatable, dimension(:,:) :: vtx_on_face
    integer, allocatable, dimension(:,:) :: vtx_on_edge
    integer, allocatable, dimension(:,:) :: edge_on_face
    real(8), allocatable, dimension(:  ) :: vtx_x
    real(8), allocatable, dimension(:  ) :: vtx_y
    real(8), allocatable, dimension(:  ) :: vtx_z
    real(8), allocatable, dimension(:  ) :: vtx_lon
    real(8), allocatable, dimension(:  ) :: vtx_lat
  contains
    procedure :: init     => icosahedron_init
    procedure :: clear    => icosahedron_clear
    procedure :: topology => icosahedron_topology
    procedure :: output   => icosahedron_output
    final icosahedron_final
  end type icosahedron_type

contains

  subroutine icosahedron_init(this, nsub)

    class(icosahedron_type), intent(inout) :: this
    integer, intent(in) :: nsub

    real(8) a, b, c
    integer i, j, n, nf, ne, isub, ivtx
    logical, allocatable :: divided(:)

    call log_notice('Initialize icosahedron with ' // to_str(nsub) // ' subdivisions.')

    call this%clear()

    this%nsub = nsub

    n = 2**nsub
    this%nv = 10 * n**2 + 2
    this%nf = 20 * n**2
    this%ne = 30 * n**2

    allocate(this%vtx_on_face (3,this%nf))
    allocate(this%vtx_on_edge (2,this%ne))
    allocate(this%edge_on_face(3,this%nf))

    allocate(this%vtx_x  (  this%nv))
    allocate(this%vtx_y  (  this%nv))
    allocate(this%vtx_z  (  this%nv))
    allocate(this%vtx_lon(  this%nv))
    allocate(this%vtx_lat(  this%nv))

    a = 1.0d0
    b = 0.0d0
    c = (1.0d0 + sqrt(5.0d0)) * 0.5d0

    this%vtx_x(1:12) = [-a,  a, -a,  a,  b,  b,  b,  b,  c,  c, -c, -c]
    this%vtx_y(1:12) = [ c,  c, -c, -c, -a,  a, -a,  a,  b,  b,  b,  b]
    this%vtx_z(1:12) = [ b,  b,  b,  b,  c,  c, -c, -c, -a,  a, -a,  a]

    do i = 1, 12
      a = sqrt(this%vtx_x(i)**2 + this%vtx_y(i)**2 + this%vtx_z(i)**2)
      this%vtx_x(i) = this%vtx_x(i) / a
      this%vtx_y(i) = this%vtx_y(i) / a
      this%vtx_z(i) = this%vtx_z(i) / a
      call inverse_cartesian_transform(this%vtx_lon(i), this%vtx_lat(i), this%vtx_x(i), this%vtx_y(i), this%vtx_z(i))
    end do

    call this%topology(0)

    ivtx = 13
    do isub = 1, nsub
      n  = 2**(isub - 1)
      nf = 20 * n**2
      ne = 30 * n**2
      allocate(divided(ne)); divided = .false.
      do i = 1, nf
        do j = 1, 3
          if (.not. divided(this%edge_on_face(j,i))) then
            divided(this%edge_on_face(j,i)) = .true.
            call middle_point(this%vtx_lon(this%vtx_on_edge(1,this%edge_on_face(j,i))), &
                              this%vtx_lat(this%vtx_on_edge(1,this%edge_on_face(j,i))), &
                              this%vtx_lon(this%vtx_on_edge(2,this%edge_on_face(j,i))), &
                              this%vtx_lat(this%vtx_on_edge(2,this%edge_on_face(j,i))), &
                              this%vtx_lon(ivtx), this%vtx_lat(ivtx))
            call cartesian_transform(this%vtx_lon(ivtx), this%vtx_lat(ivtx), &
                                     this%vtx_x(ivtx), this%vtx_y(ivtx), this%vtx_z(ivtx))
            ivtx = ivtx + 1
          end if
        end do
      end do
      deallocate(divided)
      if (isub /= nsub) call this%topology(isub)
    end do

  end subroutine icosahedron_init

  subroutine icosahedron_clear(this)

    class(icosahedron_type), intent(inout) :: this

    if (allocated(this%vtx_on_face )) deallocate(this%vtx_on_face )
    if (allocated(this%vtx_on_edge )) deallocate(this%vtx_on_edge )
    if (allocated(this%edge_on_face)) deallocate(this%edge_on_face)
    if (allocated(this%vtx_x  )) deallocate(this%vtx_x  )
    if (allocated(this%vtx_y  )) deallocate(this%vtx_y  )
    if (allocated(this%vtx_z  )) deallocate(this%vtx_z  )
    if (allocated(this%vtx_lon)) deallocate(this%vtx_lon)
    if (allocated(this%vtx_lat)) deallocate(this%vtx_lat)

  end subroutine icosahedron_clear

  subroutine icosahedron_topology(this, isub)

    use delaunay_voronoi_mod

    class(icosahedron_type), intent(inout) :: this
    integer, intent(in) :: isub

    integer i, j, n, nv, nf, ne
    type(linked_list_iterator_type) it
    type(delaunay_vertex_type), pointer :: DVT
    type(delaunay_edge_type), pointer :: DE

    n  = 2**isub
    nv = 10 * n**2 + 2
    nf = 20 * n**2
    ne = 30 * n**2
    call delaunay_voronoi_init(nv, lon=this%vtx_lon(1:nv), lat=this%vtx_lat(1:nv))
    call delaunay_triangulation(init_three_vertices=[11, 5, 2], mute=.true.)

    it = linked_list_iterator(global_DTs)
    do while (.not. it%ended())
      select type (DT => it%value)
      type is (delaunay_triangle_type)
        i = DT%id
        do j = 1, 3
          DVT => get_DVT(DT%DVTs%value_at(j))
          this%vtx_on_face(j,i) = DVT%id
          DE => get_DE(DT%DEs%value_at(j))
          this%edge_on_face(j,i) = DE%id
        end do
      end select
      call it%next()
    end do
    do i = 1, global_DEs%size
      DE => get_DE(global_DEs%value_at(i))
      this%vtx_on_edge(1,i) = DE%DVT1%id
      this%vtx_on_edge(2,i) = DE%DVT2%id
    end do

  end subroutine icosahedron_topology

  subroutine icosahedron_output(this, file_path)

    use fiona

    class(icosahedron_type), intent(inout) :: this
    character(*), intent(in) :: file_path

    call fiona_create_dataset('ico', desc='Icosahedron', file_path=file_path)
    call fiona_add_dim('ico', 'vtx', size=this%nv)
    call fiona_add_var('ico', 'vtx_lon', long_name='Vertex longitude', units='degrees_east' , dim_names=['vtx'], dtype='r8')
    call fiona_add_var('ico', 'vtx_lat', long_name='Vertex latitude' , units='degrees_north', dim_names=['vtx'], dtype='r8')
    call fiona_start_output('ico')
    call fiona_output('ico', 'vtx_lon', this%vtx_lon * deg)
    call fiona_output('ico', 'vtx_lat', this%vtx_lat * deg)
    call fiona_end_output('ico')

  end subroutine icosahedron_output

  subroutine icosahedron_final(this)

    type(icosahedron_type), intent(inout) :: this

    call this%clear()

  end subroutine icosahedron_final

end module icosahedron_mod
